rm(list=ls(all=TRUE))

library(foreign)
library(memisc)
library(xtable)
library(Hmisc)
library(maptools)
library(stringr)
library(geosphere)
library(fields)
library(MASS)
library(raster)
library(rgdal)
library(shapefiles)
library(rgeos)
library(multiwayvcov)
library(stargazer)
library(plm)
library(lmtest)
library(data.table)
library(car)
library(RColorBrewer)
#library(plyr)
#library(plotKML)
#library(ggmap)
#library(grImport)
#library(scales)
#library(readstata13)


if(Sys.info()["user"]=="IV")
  setwd("/Users/IV/")
if(Sys.info()["user"]=="benjaminbarber")
  setwd("/Users/IV/")


Combined_Data <- read.dta("/Users/IV/Dropbox/Nazi/Data/Master Data w Yearly Radio Data & ITM Final.dta")

shape <- readShapeSpatial("/Users/IV/Dropbox/nazi/nazi_data/OLoughlin Shapefiles/germany_1933.shp", 
	proj4string=crs("+proj=lcc +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs"))
shape2 <- spTransform(shape, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
shape3 <- spTransform(shape, CRS("+init=epsg:3347"))

area_test <- gArea(shape3, byid=TRUE)
area_test2 <- gArea(shape2, byid=TRUE)
cor(area_test, area_test2) # These are roughly the same...

area_test2 <- as.data.frame(area_test2)
area_test2$id <- rownames(area_test2)

soldier_dat <- Combined_Data[,c("PersID", "lon_master", "lat_master")]
length(soldier_dat[,1])

### Get rid of missing data
soldier_dat <- subset(soldier_dat, !is.na(lon_master))
length(soldier_dat[,1])

### Transform
new_data <- as.matrix(soldier_dat)
xy <- soldier_dat[,c(2,3)]

soldier_spatial <- SpatialPointsDataFrame(coords=xy, data= soldier_dat, 
	proj4string=crs("+proj=lcc +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs"))
soldier_spatial2 <- SpatialPointsDataFrame(coords=xy, data= soldier_dat, 
	proj4string=crs("+proj=longlat +ellps=WGS84 +datum=WGS84"))


Mixed <- over(soldier_spatial2, shape2)
head(Mixed)
Mixed_with_area <- merge(Mixed, area_test2, by.x="id", all.x=TRUE)


final_data <- cbind(soldier_dat, Mixed_with_area)
var_keep <- c("PersID", "population", "area_test2")
final_data <- final_data[var_keep]

Combined_Data_Pop <- merge(Combined_Data, final_data, by="PersID", all.x=TRUE)

Combined_Data_Pop$population_fix <- ifelse(Combined_Data_Pop$population<=0,NA,Combined_Data_Pop$population)
Combined_Data_Pop$log_population_fix <- log(Combined_Data_Pop$population_fix)
#Combined_Data_Pop$area_test2_fix <- ifelse(Combined_Data_Pop$population==0,NA,Combined_Data_Pop$area_test2)

#Combined_Data_Pop$pop_density_guess <- Combined_Data_Pop$population_fix* Combined_Data_Pop$area_test2_fix

Combined_Data_Pop_Export <- Combined_Data_Pop

indx2 <- which(sapply(Combined_Data_Pop_Export, is.character))
for(j in indx2) set(Combined_Data_Pop_Export, i = grep("^$|^ $", Combined_Data_Pop[[j]]), j=j, value=NA_character_)


write.dta(Combined_Data_Pop, "/Users/IV/Dropbox/Nazi/Data/Master Data w Yearly Radio Data & ITM with Pop Final.dta")

